Efficient kinetic Monte Carlo method for reaction-diffusion processes 
with spatially varying annihilation rates 
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We present an efficient Monte Carlo method to simulate reaction-diffusion processes with spatially 
varying particle annihilation or transformation rates as it occurs for instance in the context of 
motor-driven intracellular transport. Like Green's function reaction dynamics and first-passage 
time methods, our algorithm avoids small diffusive hops by propagating sufficiently distant particles 
in large hops to the boundaries of protective domains. Since for spatially varying annihilation or 
transformation rates the single particle diffusion propagator is not known analytically, we present 
an algorithm that generates efficiently either particle displacements or annihilations with the correct 
statistics, as we prove rigorously. The numerical efficiency of the algorithm is demonstrated with 
an illustrative example. 
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I. INTRODUCTION 

Kinetic Monte Carlo simulations are frequently used in 
various fields to analyze the spatio-temporal evolution of 
systems consisting of many freely diffusing particles that 
can collide, react, transform or annihilate. Spatial as 
well as stochastic aspects are important when diffusion 
is not sufficiently fast to make the system well-stirred 
and the number of reactants within diffusion range is 
small. In this case a mean-field description, for instance 
with a set of coupled reaction-diffusion equations, is in- 
appropriate. Moreover, in the limit of extreme dilution 
methods using a discretization of the underlying stochas- 
tic reaction-diffusion system, either in time [l| or in space 
[HEl: become computationally inefficient. 

The currently most efficient methods to simulate ex- 
tremely diluted reaction-diffusion systems are Green's 
function reaction dynamics Jjjp and first-passage ki- 
netic Monte Carlo methods In essence they avoid 
the small diffusion hops of a conventional random walk or 
Brownian dynamics simulation by propagating particles 
over long distances through a sequence of large displace- 
ments. The latter are generated stochastically according 
to the exactly known Green's function for a freely diffus- 
ing particle within so-called protective domains that are 
free from other particles. The typical size of these pro- 
tective domains is inversely proportional to the particle 
density and the larger these domains are (i.e. the smaller 
the particle density is) the more efficient the algorithm 
is. 

In general, during the free diffusion the particle can 
also be annihilated or transformed with a rate k into a 
different species, in which case the Green's function is 
still exactly know. In this paper we address the question 
how to propagate the particles when the annihilation rate 
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varies in space and time, denoted as fc(r, t). This prob- 
lem arises for instance in the context of motor-driven 
intracellular transport, where particles (or cargos) can 
in addition to diffusion and reaction also attach to a cy- 
toskeleton filament and move ballistically with a constant 
speed in the direction of the filament. A continuum de- 
scription of the diffusive and ballistic modes of motion 
[H, [l(| involves the filament density p(r, t) which deter- 
mines the local rate with which freely diffusing particles 
make a transition into the ballistic state. In a typical 
cell the filament density is spatially inhomogeneous and 
thus has to be taken into account during the propagation 
of particles on large scales. Analogous examples arise in 
systems in which the annihilation of particles depends 
on a spatially inhomogeneous concentration field of an 
abundant reaction partner (i.e. whose density is much 
larger such that a continuum description is appropriate 
for it). 

Green's function reaction dynamics and first-passage 
time Monte Carlo methods reduce the simulation of 
a many-particle reaction-diffusion system to individual 
particles that diffuse freely as long as other particles are 
sufficiently distant (i.e. outside the interaction range), 
and perform a reaction event once a particle pair reaches 
a minimum distance. Algorithmically one can ensure free 
diffusion for instance by estimating the maximum diffu- 
sion distance [1, Q until a reaction is scheduled or by 
the definition of protective domains for each particle p- 
Q depending on the actual arrangement of neighboring 
particles. In both cases one then utilizes the free diffu- 
sion propagator within predefined domains to generate 
stochastically a time when either the maximum distance 
is undcrrun or a protective domain boundary is reached. 
For free diffusion this is achieved using the analytically 
known Green's function, but for free diffusion with spa- 
tially varying annihilation rates this propagator is unfor- 
tunately not analytically available. 

Thus in this paper we consider a freely diffusing single 
particle in an arbitrary domain G € R n that can be an- 
nihilated with a time and space dependent rate k(r, t). 
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In general, annihilation means a transition into a differ- 
ent species that is not considered in the present reduced 
setup. For a particle initially at time to located at ro G G 
this diffusion-annihilation process is described by the fol- 
lowing diffusion-annihilation equation 



dP(r,t\r ,t ) 

at 



DAP(r, t\r , t ) - k (r, t) P(r, t|r 0> t ) , 



(1) 



where P(r, t\ro,to) is the probability density to find the 
particle at time t at r G G. For arbitrary k(r, t) and arbi- 
trary G there is no analytic solution of Eq. (JTJ available. 
In principle this equation can be solved numerically, but 
in the context of a general reaction-diffusion system (in- 
volving many particles and several particle species) using 
for instance the first-passage Monte Carlo method this is 
unfeasible: Here one needs for each particle hop the whole 
first-passage time distribution for a particle to reach the 
protective domain boundary dG, which is computation- 
ally too demanding to be carried out in the innermost 
loop of the algorithm. 

Therefore we present in this paper an algorithm that 
samples times t > to and positions r for arbitrary anni- 
hilation rates k(r,t) and arbitrary domains for which a 
particle diffusing according to Eq. (JTJ) either a) reaches 
the boundary for the first time (r € dG) or b) is anni- 
hilated (r G G). In addition, a slightly modified version 
of the algorithm generates the whole probability density 
P(r, i|ro, to) within G, meaning it solves Eq. (JTJ) stochas- 
tically. 

The paper is organized as follows: Section [TT] defines 
all probability densities and flows used throughout this 
paper. Based on the ideas of 043, section|nT]presents an 
algorithm for the sampling of (r, t) on arbitrary domains 
G in the case of a spatially homogeneous but tempo- 
rally varying annihilation rate k(r,t) = k(t). Section HVl 
generalizes this method to spatially inhomogcneous rates 
k (r, t), proves its correctness and discusses its efficiency. 
Finally section [V] shows an application example of this 
method. 



II. DEFINITIONS 

In this section the probability densities and flows used 
later on are defined. Let P(r, t\ro, to) be the probabil- 
ity density solving the diffusion-annihilation equation (JTJ) 
within the domain G with boundary dG, possibly partly 
absorbing, partly reflecting. The particle annihilation 
generates a probability flow / a (r, t\ro, to) out of the sys- 
tem given by 



f a (r,t\v 0l to) = k(r,t) ■ P(r,t\r ,t ) . 



(2) 



The probability flow /b(r, t\ro, to) at the absorbing 
parts of the boundary at time t at position r € dG is 
given by 



f b (r,t\r ,t ) = -D VP(r,t|r 0j to) ■ n r 



(3) 



where n r denotes the outward pointing unity vector 
perpendicular to the boundary dG at r. Consequently 
P(r, t\vo,to) is not normalized for t > to- The corre- 
sponding probability density p e (t\ro,to) for an annihila- 
tion or absorption event is given by 



Pe(t\r ,t ) 



d 
~dt 



I dr P(r,i|r ,io) 
Jg 



a(t|r ,to)+j8(f|r 0j to) 



(4) 



with a(t\r ,t Q ) = / dr f a (r, t\r , t ) 
Jg 



and /3(t\r ,t ) 



i)G 



dF f b {r,t\r ,t ) 



where dF denotes the surface element at position r G dG. 
Hence, the task is to sample the pairs (r, t) in statisti- 
cal agreement to / a (r, £|r , to) and fb(r,t\r ,to), i.c the 
statistic of t will be according to p e . 

In the following we also need the probability distribu- 
tion of a freely diffusing particle Pq (r, i|ro,£o) without 
annihilation, which obeys 



dPp (r,t\r ,t 

dt 



D AP D (r,t\r ,t Q ) 



(5) 



The probability density for being absorbed at the 
boundary for a purely diffusing particle at time t is given 
by 

d f 

p b D {t\r ,t ) = - — J^dr P D (r,t\v ,to) (6) 

and the probability density of the absorbing position r G 
dG under the condition that the absorption takes place 
at time t is given by 



pf (r|i,r ,io) = 



VPD(r,t|r ,to)-n r 
J dG dF VPd (r,t|r ,to) • n r 



(7) 



Using the Gauss's theorem and Eq. ([5]) in the denomi- 
nator, one obtains: 



pf(r|t,r ,fo) = -D 



VP D (M|r ,t )-n r 
Pb(t\r ,to) 



(8) 



For spatially homogeneous annihilation rates fc(r, t) = 
k(t) the annihilation process decouples from all spatial 
variables, i.e. the solution of Eq. (JlJ can be written as 



P(r,f|r 0l io)=e- / 'o fe (*') dt '.p D (r,t|ro,i ) 



(9) 



Hence, the probability density of being annihilated at 
time t under the condition of not being absorbed at the 
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boundary before for a spatially homogeneous rate k(t) is 
given by 



- f* k(t')dt' 



(10) 



and the probability density of the annihilation position 
r £ G under the condition that the particle is annihilated 
at time t is given by 



p^(r|t,r ,t ) = 



P D (r,t\r ,t ) 
f G dr> P D (r>,t\r ,t ) ' 



(11) 



which is equal to the probability density of a purely dif- 
fusing particle under the condition of not being absorbed. 



III. HOMOGENEOUS ANNIHILATION RATE 

In this section we present an algorithm that samples 
times t > to and positions r for homogeneous annihi- 
lation rates k(r,t) = fc(t) and arbitrary domains G for 
which a particle diffusing according to Eq. ([TJ either a) 
reaches the boundary for the first time (r £ dG) or b) is 
annihilated (re G). 

Assuming the ability to generate random numbers ac- 
cording to all the densities, which were defined in the 
previous section, a correct way of sampling (r, t) is shown 
in Algorithm [TJ 

Algorithm 1 KMC 1 

Input: r , t 
Output: r, t 

t a random number according to Pa{-\to) 
tb <— random number according to (-|ro,to) 
t •(— min(t a , tb) 
if (t a < h) then 

r random position according to p®(-\t a , ro, to) 
else 

r random position at the boundary dG 
according to pf(-\tb, r , to) 

end if 
return (r, t) 

The probability density A(r, i) that the algorithm pro- 
duces an annihilation at time t at position r is then given 
by 



A(r,t) = (g(t\ta)U dt b pf(t b \r ,to)j ^(r|t,r ,t ) 

= p%(t\t )P D (r,t\r ,t ) 

= fa (r,t|r ,t ) • 

The probability density B(r,t) that the algorithm de- 
livers an absorption at time t at position r <G dG is given 
by 

B(r,t) = (J dt a p°(t a \t )\ p£(t\r ,t )pf(r\t,T ,t ) 
= fb(r,t\r ,t ) ■ 



Consequently, the statistic of random pairs (r, t) gen- 
erated in this way coincides with f a (r, t|r , to) and 
/b(r, t|ro, to) and is therefore correct. 

One problem remains: As there are no analytic solu- 
tions for Eq. ([3]) available for arbitrary domains G, it 
is not possible to sample the quantities p®, p® and p® 
directly for arbitrary domains G and arbitrary boundary 
conditions. Only a direct sampling of p® is possible, as 
there is always an analytic expression for the correspond- 
ing distribution function available: 



F^{t\t ) = l-e-^ k(t ' )dt ' 



(12) 



That means, Algorithm [TJ is useful only in some special 
geometries G. 

Nevertheless, it is possible to use these special cases 
to sample the random pair (r, t) for arbitrary domains. 
Two different methods will be shown now. 



A. Subset method 

In a kinetic Monte Carlo method for the 

simulation of reaction-diffusion processes of many-body 
systems is presented. It is based on the fact that there 
are analytic solutions of Eq. ([5]) for some simple domains 
G' and boundary conditions. The appendix shows a 
list with some of these domains in one, two and three 
dimensions and derives expressions for distribution 
functions, which are necessary for the usage of the 
inversion method (Tlj . 

If G' denotes a subset of G with ro £ G", the shape 
of G \ G' will not matter for the particle, as long as the 
particle has not left G' for the first time. Hence, we can 
restrict the description of the particle's motion to G' un- 
til it leaves G' for the first time. Mathematically we are 
dealing with a first-passagc-problem in G' . Its solution is 
given by Eq. ([5]) on G' according to absorbing boundary 
conditions at the interior of G and the boundary condi- 
tions of G at common boundaries of G and G' (as far as 
they exist). 

Assuming that we arc able to sample all occurring 
densities, a random pair (r,t) for G' can be generated, 
as shown in Algorithm [JJ If annihilation takes place 
(t Q < tfc), the particle is annihilated before it leaves G' 
and therefore it is not influenced by the restriction to 
G' . If the particle reaches the boundary of Q (t < t a ), 
two possibilities have to be distinguished: For r £ dG it 
reaches an absorbing boundary of G and the algorithm 
will stop. Otherwise, the particle continues its diffusive 
motion under the condition of having been at position 
ro = r at time to = t. As it is always possible (see ap- 
pendix) for an arbitrary ro to find a subset of G where 
there are all needed analytic expressions available, we 
can go on this way until the particle is annihilated or 
absorbed at the boundary of G. The pseudo-code of this 
is shown in Algorithm [21 
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Algorithm 2 KMC 2 



Input: r , to 
Output: r, t 

t «- to 
r «- r 
repeat 

choose a suitable domain G' with r 6 G' 
t a <— random number according to pf (-|t) on G' 
tb <— random number according to pf (-|r, t) on G' 
if (t a < tb) then 

r <— rand, position according to (-\t a ,r,t) on G' 
else 

r <— random position at the boundary 9G' 
according to pf (-\tb,r,t) 

end if 

t <- min(t a , t b ) 
until ( r G <9G or t a < t b ) 
return (r, t) 





FIG. 1: Illustration of the method: Absorbing boundaries are 
shown in green, reflecting ones in red. In all situations a)-d) 
the case t a > tb is sketched, otherwise the algorithm would 
stop earlier. Depending on the position of the particle the 
shape of the chosen domain G' varies (rectangles and circles) . 



A sketch of the method in a case where the particle is 
absorbed at the boundary is shown in Fig. [TJ 

For a given domain G, the efficiency of the method 
depends on the choices of G". Ideally, one chooses G' 
from the list of possibilities in a way that maximizes the 
expectation value of t b . However, it also takes more time 
to look for this special subset and eventually calculate 
the random numbers for this situation. In the cases of a 
particle in the middle of a circle (sphere) or in the middle 
of a square (cube) the random numbers can be generated 
very fast. Hence, in some situations it might be better 
to use smaller domains G' than in principle possible. 

Up to this point, there is no approximation involved, 
but depending on the shape of dG a problem occurs: If 
the particle approaches an absorbing part of the bound- 
ary of G, it will always automatically approach an ab- 
sorbing part of the boundary of the chosen G', too, as 
G' C G. In consequence, the expectation value of tb will 
decrease, the stopping condition t a < tb becomes more 
and more unlikely and the time incrementations in t will 
become smaller and smaller, if the condition r € dG is 
not fulfilled. But r £ dG can only be true, if the intersec- 
tion of dG and dG' contains more than just single points. 



The same problem occurs for reflecting boundaries of G. 
As the choice of G' is limited, we sometimes have to ap- 
proximate dG by a polygon in order to avoid a critical 
slowing down of the algorithm. However, it is important 
to mention that we can always choose the accuracy of the 
approximation by the choice of the polygon. 



B. Maximum distance method 

Depending on the shape of dG close to the position 
of the particle, there is sometimes a better way of prop- 
agating the particle than it is shown in the subsection 
above. [J, [f| introduced this idea for the particle's short 
time behavior in the context of particle-particle interac- 
tion, but it can be modified for a usage in our context. 
It is based on the assumption that there is a maximum 
distance Ar, which the particle does not reach within 
a time At. Hence, within this time interval At, only 
the intersection of G with a neighborhood of radius Ar 
matters. Of course, this assumption is an approximation 
since there is a non-vanishing probability that the par- 
ticle leaves this neighborhood within At. However, it is 
possible to control the accuracy by the definition of Ar 
via a parameter 7. We define: 



Ar = 7VD At 



(13) 



As 7 increases, it is more and more unlikely for the 
particle to violate the assumption. More precisely it is 
even possible to give an upper boundary for failing the 
assumption by studying the first passage-process to the 
boundary of a particle that starts in the middle of a circle 
(2d) or a sphere (3d) with radius Ar and calculating the 
probability ^(7) for not having reached the boundary 
within At. 

In 2 dimensions we obtain: 



w 2 d(l) = 2 5Z 



1 



1 



^— ' a„ Ji(a„) 

n— 1 v 1 



e > 



(14) 



where a n , n € N are the roots of the Bessel function 
Jo (see appendix and be aware of the slightly different 
notation). In the following tabular the corresponding 
values are calculated for some 7. 



7 


2 


3 


4 


6 


7 


9 


1 - W 2 d(l) 


0.623 


0.193 


0.0347 


2.41e-4 


9.39e-06 


3.90e-9 



In 3 dimensions we get (see appendix): 

00 

i(7) = 2^(-l)»+ 1 e -(^) 2 



W 3 d{ 



(15) 



In the following tabular the corresponding values are 
calculated for some 7. 
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FIG. 2: For t < Ar 2 /^ 2 D, the particle is assumed not to 
cross the black circle line. So the solution of an infinite sector 
with reflecting boundaries can be used. 



new (r, t). 

For all t > to we define the spatially homogeneous but 
time dependent upper bound for the annihilation rates 



k m (t) = max reG {fc(r,t)} . (16) 

The density p® with the rate k m (t) is denoted by p m 
in the following: 



7 


2 


3 


4 


6 


7 


9 


1 - w 3d (l) 


0.830 


0.357 


0.0827 


8.36e-4 


3.78c-05 


1.63e-8 



Pm{t\ta) = - — 



- f! k m (t')dt' 



(17) 



Consequently, for a choice of 7 in the range of 7 — 9 
one is on the safe side for all practical purposes, where 
also other numerical error sources (quality of the random 
number generator, rounding errors) come into play. 

This gives the possibility to use analytic solutions of 
Eq. ([5]) of domains which coincide with G only in the 
neighborhood of Ar. The example of Fig. [2] shows the 
left part of the domain from Fig. [TJ Choosing Ar in 
the shown way, the analytically known solution of an 
infinite sector (see appendix) can be used, as long as 
t < Ar 2 / 7 2 D. 

Hence, if the particle is neither annihilated nor ab- 
sorbed within At, the particle will stay diffusive and a 
random pair (r, At) must be created for the particle's 
new position. In order to avoid repetitions, we skip the 
pseudo-code details here, as they will be shown in the 
next section in a more general case. 

In most situations it is much better to use the subset 
method as its time- increments are generally much larger. 
But in a situation like the one sketched in Fig. [2j the 
particle is very close to the reflecting boundaries and no 
suitable large domain G' is available. In consequence, 
G' would be very small, leading to a very small tj, on 
average. 



IV. INHOMOGENEOUS ANNIHILATION RATE 

The last section showed how to find a solution for an ar- 
bitrary domain G by solving the problem in several steps 
in smaller domains G". Hence, without loss of generality, 
we now assume the ability to sample random numbers 
according to pf, p® and p® directly. 

If the annihilation rate becomes inhomogeneous, Eq. 
© is not a solution of ([T]) anymore. The annihilation- 
time is now correlated to the complete path of the par- 
ticle, thus the method presented in the previous section 
will not work. In this section we present a way to over- 
come this problem for arbitrary rates k(r, t) without any 
additional approximations. The following method starts 
with the pair (ro,to) and generates a series of random 
pairs (rj,£j). The last pair of this series will become the 



We sample a candidate pair (ri,ti) as shown in Al- 
gorithm [T] For ri G dG the particle is absorbed at the 
boundary, i.e. the first candidate is accepted. Otherwise 
we compare the ratio k(ji,ti)/k m (ti) to a uniformly dis- 
tributed random number a; in [0,1]. If k(ri, ti)/k m {t\) > 
x, the particle is annihilated, i.e. the first candidate 
is also accepted, else we generate a new candidate pair 
(jC2, £2) under the condition of having been at position ri 
at time t\. This can be continued until the particle is 
absorbed at the boundary of G or annihilated. 

The algorithm can also be used to sample the com- 
plete probability density P(r, i|ro, to): If no candidate 
is accepted until an arbitrarily chosen time t max = t is 
reached, the algorithm returns a random position of the 
still diffusive particle, i.e a pair (r,t max ) whose statistics 
is given by P(r, t max \r , t ). Also in case one wants to 
use the maximum distance method, the time t max has to 
be chosen appropriately If a break at t max is not wanted, 
one simply sets t max — 00. A pseudo-code description is 
shown in Algorithm [3] 

Algorithm 3 KMC 3 

Input: r , to, tmax, k m (t) 

Output: r, t 

t <- to 

r <- r 

repeat 

t a <— random number according to p m (-\t) 
tt random number according to (-|r, t) 
if (tmax < min(t a ,t b )) then 

r <— random position according to p°(-|t ma x, r, t) 

t ^ tmax 

else 

if (t a < tb) then 

r <— random position according to p^(-|t a , r, t) 
else 

r random position at the boundary dG 
according to pf (-|t(,, r, t) 

end if 

t min(t a , t b ) 
end if 

until > ran[0, 1]) or (t a > h) or (t = t max )) 

return (r, t) 
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A. Proof of correctness of Algorithm [31 



Starting with i = 0, we compute Wi, Ai, Bf. 



The basic mechanism by which the algorithm handles 
a spatially varying annihilation rate fc(r) is to generate 
trial annihilation positions using the propagator for a 
spatially constant (but maximal) annihilation rate k m . 
The annihilation is then accepted with the local probabil- 
ity k(r)/k m for. At first sight it appears counter-intuitive 
that this local procedure actually gives the correct statis- 
tics, since the probability to propagate a particle from ro 
to r depends on the complete annihilation rate landscape 
in between and around. Why is it sufficient to probe fc(r) 
locally at one or a few positions generated by the algo- 
rithm? 

Before we answer this question rigorously by proving 
that it is indeed sufficient, we give an intuitive argument 
why one might expect the procedure to be correct: The 
stronger the spatial variation of fc(r) is in G the larger the 
maximum rate k m will be. A large constant annihilation 
rate k m gives rise to a particle propagator that forbids 
large hops, which implies that the algorithm will produce 
many small intermediate hops and after each hop evalu- 
ates fc(r). In this way the algorithm explores stochasti- 
cally the annihilation landscape on finer or coarser length 
scales depending on the variation of k(r). If for instance 
fc(r) = everywhere in G except in a small restricted 
region, where it is fc(r) = k\ > 0, thus k m = k\. Then 
the algorithm will explore the complete region G with 
a hop size that is characteristic for the restricted region 
with the non-vanishing annihilation rate. In the end this 
yields the correct statistics for the whole region, which 
we will prove now. 

We will prove that the statistic of the output pairs 
(r,t) satisfy the probability flows f a and /& for t < i max . 
Then the case t = t max (particle is still diffusive at time 
imax) occurs with the correct probability, too. We also 
prove that the statistic of output pairs (r, i max ) coincides 
with P(r,i max |r ,t ). 

The algorithm will stop after a (unknown) number i+1 
(i € N) of loop-runs (see Algorithm^. The probability 
density for being annihilated after i + 1 loop- runs at time 
t at the position r is denoted by Ai(r,t\r ,t ). Analo- 
gously the probability density for being absorbed at the 
boundary after i + 1 loop-runs at time t at the position 
r is denoted by Bj(r, £|ro, fo)- The probability density 
for stopping after i+1 loop-runs, still being in the dif- 
fusive state at t maK and being located at r is denoted 
by Wj(r, i max |r , £o)- As the number of loop-runs is a 
disjoint decomposition, we can sum over i to obtain the 
total densities for the corresponding events: 



A(r,i|r ,i ) 
B(r,f|r 0) io) 
W{r,t max \rQ,to) 



J2Mr,t\r ,t ) , (18) 

i=0 
oo 

J2Bi(r,t\r ,t ) , (19) 

oo 

X^MmaxIro^o) ■ (20) 



= 0: 

Ao: For this event t a has to be smaller than tb, which 
delivers the second factor in the following product. The 
third factor belongs to the choice of the position and the 
last one arises from the exit-condition of the algorithm's 
loop: 



A (r,t\r ,t ) = p m {t\t ) 



dt b Pb(tb\r ,t ) 



P D (r,t\r ,t ) k(r,t) 



f G dr'P D (T>,t\r ,t ) k m (t) 
= k(r,t)e- I > t o k ^ dt 'p D (r,t\r Q ,t ) . (21) 

Bq: An analogous procedure delivers 

B (v, t\r , t Q ) = pf (t\r , to)e~ %> km[t ' )dt ' 
V(Pg(r,t|r ,t )) n r 
' J dG dF V(P D (r,t\r ,t ))-n r 

= -De-^ k ^ dt 'vP D {v,t\v 0l t )-n r . (22) 



Wo'. If the particle reaches the time i max in the first 
loop-run, t a and tb have to be larger than £ max . Thus 
Wo is the product of these two independent probabili- 
ties with the spatial density p^(r, t max |ro, to) : 



W (r,t 

max |r ,io) = 
dt b Pb(tb\r ,to) 



dt a Prn{ta\ta) 
max 



e J *o 



J G dr'PD(r',i max |ro,£o) 

fir* k ^(t')dt' : 



'max |r ,t ) ■ (23) 



1: 



As the algorithm will pass the loop twice here, we have 
to sum/integrate over all weighted pairs (ri,£i), which 
will be achieved in the first loop-run. Since the algo- 
rithm will only continue with a new loop if t a is smaller 
than tb, for the first loop the factors and integrals look 
the same for all cases. The factors of the final loop can 
be taken from the individual factors of i = with the 
starting position ri and the time t\ instead of ro and 
to- Defining the probability that the algorithm denies 
a candidate pair (r, t) 



q(v,t) =11- 



(24) 



we get: 



i=0 
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Ai(r,i|r ,to) = / dt al p m (t al \t Q ) [ dt bl p b (t bl \r ,t ) [ dr 1 ^j^rr^T^Z , q{ri,t al ) A (r 1 t\r 1 ,t al ) 

Jt Jt al JG J G ">T "d{? ,ta,l\To,to) 

= k{v,t)e- S ^ km[t ' )dt ' [ dt a i f dr 1 P D {r 1 ,t al \r ,t )k m {t al )q{r 1 ,t al )P D {v,t\T 1 ,t al ) (25) 
J t Jg 

Bi(r,i|ro,io) = / dt al p m (t al \t ) [ dt bl p b (t bl \r ,t ) [ dr x f ^j^prpr^ 7 g(ri,t a i) Po(r,t|n,t a i) 

Jt Jtal 



-D e J *o 



-ft k m (t')dt' 



G f G dr'P D (r',t al \r ,t ) 

/ dt al / driPi3(ri,i a i|ro,to)fc m (t i)g'(ri,i i)VPD(r,i|ri,i a i) • n r (26) 
7tn Jg 



VFi(r,t max |r ,t ) (27) 

d*al Pm(tal|*o) I 00 dt bl P ?(t bl \r ,t ) / f P f^ a ^° M \ n{v u t al )W {v,t 

max r l; *al) 

/to Jt al Jg J G dr'P D (r',t a i|r ,t ) 

= e~^'o fe ™(*) dt / / P D (ri,tai\ro,to)k m (tai)q(ri,tai)PD(r,t max \ri,tai) (28) 

For i > 1 we introduce the definitions 



(29) 



=0..i) ■f > D( r i*l r i)*i) ' 
t <ti<...<ti<i 

with Qi({( r i^0}i=o..i) = n i5 D( r '' i 'l r i-i'*'-i)-9(ri,ii) ■ (30) 

2=1 

Finally, defining 

fto (r, t|r , to) = Pd (r, t|r , to) , (31) 

one inductively gets for i > 0: 

A 4 (r,t|r ,t ) = k(r,t)e- f * t ° k ' n ( t ' )dt 'h i (T,t\ro,t ) , (32) 
BiMro.to) = -iJe-^^'^V/hCr.tlro.tb)-!!, , (33) 
Wi(r,t max \r ,t ) = e- S ^ kmit ' )dt 'hi(r,t max \TQ,t Q ) . (34) 

Hence, the total probability densities can be written as 

A(r,t\r ,t ) = k(r,t) P(r,t\r ,t ) , (35) 

B(r,i|r ,i ) = —D VP (r, t\r ,t ) ■ n r , (36) 

W(r,t max \r ,t ) = P(r,t max \r ,t ) , (37) 

with 

oo 

P(T,t\r ,t ) = e- f 'o km(t ' )dt 'j^hi(r,t\T ,to) . (38) 

i=0 

Comparing this with the definitions of / and it remains to show that P (r, t\ro, to) = P (r, i|ro, to), i.e. P (r, t|ro, to) 
has to satisfy Eq. ([T]) with the initial condition P(r, i l r o j to) = <K r — r o)- As all ft,,; with i > 1 vanish for t = t , the 
initial condition is simply fulfilled by the definition of Iiq. For the time-derivative of hi one inductively gets for i > 1: 
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hi(r,t\r ,t ) = k m (t)q(r,t) ■ hi-i(r,t\r ,t ) 



-J ljl dt n km ^J Gn (n/ rfe ) ^^'^1=0.4) PDir^ti) 



(39) 



t <t 1 <...<t i <t 

Hence, the time-derivative of P (r, f |ro, to) satisfies 



dP(r,tJr ,t ) = p ^ ^^^^ + e _ Sia km(t > )dt > I p D ^ t|ro> + g ^ ^ 



(40) 



f OO 

00 « / i \ C ( ^ 

/ n*i / n« 

i=i y=i y Jg " \fe=i 



dr k Qt ({{rhti)} l=0 J P D {r,t\ri,U) 



t a <t 1 <...<t i <t 



(41) 



-k(r, t)P (r, t|r , to) + e ^ fcm(t ' )d *' {p c (r, t|r , t ) + 



£/ (n^j fcmfe) / G „ (n**) Qj 



({(^OW^M^) 



t <ti<. ..<*;<* 

Thus, using Eq. ([5]), it follows: 
0P(r,t|r o ,<o) 



(42) 



9t 



fc(r,t)P(r,t|r ,t )+e' 



- f» tmax k m (t')dt' 



OA (r,t|r ,t ) 



. i=0 



= -fc(r, t)P (r, t|r ,*o) + D AP(r, t|r , t ) , 
which is exactly Eq. (jTJ) . Hence, the correctness of the Algorithm [3] is proven. 



(43) 



V. EXAMPLE 



This section presents a two-dimensional application ex- 
ample of the algorithm. It is designed to demonstrate 
how the algorithm handles a situation in which its cor- 
rectness is most counter-intuitive: We choose the anni- 
hilation rate to be non- vanishing just inside a restricted 
region, a circle, where it oscillates in time and varies spa- 
tially. For a chosen test-setup, we compare its results 
with the solution of a commercial FEM (finite element 
method) routine. 

At to = the diffusing particle (D = 1) is located at 
the position ro = (0; 5) within a rectangle of size 10 x 5. 
The right boundary is chosen to be absorbing, all other 
boundaries are reflecting. A strongly anisotropic time 
dependent annihilation rate fc(r,t) is chosen to be 



fc( r)i ) = / 3 l cos3 (l)|-( c2 -|l r - z l| 2 )' H r - Z ll< c ,(44) 
I , ||r-z|| > c 

with z = (5; 1.25) and c = 1.25. Fig. [3]presents a sketch 
of the described setup. 

On the one hand the problem has been solved numer- 
ically by applying a commercial FEM solver with a very 
fine triangulation (> 600000 elements) to Eq. (fTJ). In 
the following this solution is denoted by Pp(r, t). On 
the other hand the Monte-Carlo algorithm has been ap- 
plied to the problem in 4.2 x 10 8 samples. In principle 
it is not necessary to use the subset method here, as the 
analytic solution of Eq. ([5]) is known for the rectangle 
(see appendix), from which all occurring densities can be 
sampled. Nevertheless it has been used, as it increases 
the speed of the algorithm dramatically: For all highly 
anisotropic annihilation rates the ratio k/k m in Algo- 
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FIG. 3: sketch of the simulation setup: The particle starts 
its diffusive motion at the upper left corner. It can either be 
absorbed at the right wall or annihilated within the drawn 
circle. 



rithm [3] will mostly be very small (in our case even 0). 
Hence, a lot of loop runs with just small time incrementa- 
tions will on average be needed for an event. Restricting 
the movement of the particle temporally to a subset of 
G (subset method) ensures the possibility of choosing a 
smaller k m (t), which reduces the number of loop- runs in 
Algorithm [3] dramatically. Using this in our case it takes 
around 40 minutes for 10 6 samples on a single core with 
3.4 GHz. 

Firstly, the relative frequencies for the times of an 
event and the kind of the event were counted. Smc{^) 
denotes the relative frequency of having had no event 
until time t. It has to be compared with the value of 
Sp(t) = 1 — J p e (t'\ro,to)dt', which was derived numer- 
ically from the FEM solution. AMc(t) denotes the rel- 
ative frequency of having had an annihilation event be- 
fore time t. It is compared to Ap{t) = J * a(t'\ro, to)dt' . 
BMc(t) denotes the relative frequency of being absorbed 
at the right boundary before time t. It is compared to 
Bp(t) = f Q P(t'\ro,to)dt' . A plot of these quantities is 
shown in Fig. [4] 




10 20 30 40 50 60 70 80 90 100 

t 

FIG. 4: Comparison of the relative frequencies derived from 
Algorithm [3] and the numerically integrated probabilities for 
the events of still being diffusive (S), having already been 
annihilated (A) and having already been absorbed at the right 
wall (B). 



It shows an almost perfect coincidence of all corre- 
sponding quantities. The maximum relative deviation is 
about 1% in all curves. Keeping in mind that SV(i), 
Ap(t) and Bp(t) are calculated by a numerical time- 
integration of a numerical spatial integration of a nu- 
merical solution of Eq. ([lj, these small deviations are 
explainable. More precisely, Sp(t) + Ap{t) + Bp(t) = 1 
has to hold for all times, but the numerical discrepancy 
in this sum is also about 1% at maximum. 

Secondly, we want to compare the spatial distribution 
of the particle's position from the KMC algorithm to 
Pp(r,t) for three characteristic times: t\ = 5, t<z = 10, 
ts = 20. Hence, the rectangle is divided in 2N x N 
squares s xy (x = i- jj, y = j • j^, i = 0, 1, . . . 2N- 1, j = 
0, 1, . . . N— 1) and the relative frequency h xy (t) for being 
at the square s xy is counted for t\, t 2 and £3. Technically 
this has been done by setting t max = ti (i 6 {1, 2, 3}) in 
Algorithm |3] The quotient of h xy it) and the area of a 
square element is denoted by P xy (t). This density con- 
verges to the solution of Eq. ([!} in the limits of increasing 
sample numbers and TV — > 00. The upper panel in Fig. m 
shows the density P xy (t) for the chosen times in a 3d plot 
for N — 50. In order to illustrate the influence of the an- 
nihilation within the circle, the projection on the bottom 
shows isolines by discretising the density into intervals. 

• t\ = 5: The probability density of the particle is still 
centered around the starting position in the upper left 
corner. Nevertheless its shape is already influenced by 
the annihilation within the circle. 

• t 2 = 10: At time t 2 there has been almost no annihila- 
tion for a short period (slope of the red line in Fig. @|. 
Hence, diffusion almost equilibrated the density gradi- 
ent in y-direction, generated by the annihilation within 
the time-interval [5; 8]. 

• £3 = 20: At time £3 relatively strong annihilation takes 



place, which even 
within the circle. 



leads to a local minimum of P. 



In order to quantify the local differences between the 
KMC result and the FEM result, we choose squares of 
size 0.2 x 0.2 (N ~ 25). A measure for the spatially 
resolved relative deviation is 



h xy (t)-f: +0 - 2 d X >r^dyip F (T>,t) 



y+0.2 



f: +0 - 2 d X >n + »Uy>P F (r>,t) 



The lower panel of Fig. [5] shows A xy (t). For all times 
the deviations are small and in the range of the numerical 
expectation: For 4.2 x 10 s samples and (2A) 2 = 10 4 
plaquettes one expects ca. 10 4 samples per plaquette 
and thus statistical fluctuation of the order of 10 2 , i.e. 
relative fluctuation in the 1 percent range, which is what 
the lower panel of Fig. [5] confirms. 

On the right side of the simulation rectangle, where 
the absorbing boundary is located, the statistical error is 
larger for small times, since the density is still centered 
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FIG. 5: Top: The probability density P xy {t) — P(r,t\ro,to) for the two-dimensional geometry depicted in Fig. [3] and the 
annihilation rate given by Eq. (|44|) generated by the intermediate positions of Algorithm [3] for the times ti = 5, ti = 10 and 
t$ — 20. The position of the annihilation zone with oscillating strength is indicated by the full circle. Bottom: Relative 
difference A xy (t) between the Monte Carlo result P xy (t) and the numerical solution of the corresponding annihilation-diffusion 
Eq. HJ for the times h = 5 , t 2 = 10 and t 3 = 20. 



around the starting point in the upper left corner, giving 
a region (x > 8) with very small values of h xy (t). But 
the aim of the Algorithm [3] is not the stochastic solution 
of Eq. ([T]), for which finite element methods are suit- 
able. The aim of the algorithm is to sample correctly next 
events (annihilations or first-passages of boundaries) ac- 
cording to Eq. ([1]), which can not be handled by a FEM 
routine. The example demonstrates, that this is possible, 
even in cases of highly anisotropic and time dependent 
annihilation rates. 



VI. DISCUSSION 

Wc have presented an algorithm that samples correctly 
the probability distribution of a diffusing particle with 
a space dependent annihilation or transformation rate 
k(r) for arbitrary domains. Together with first-passage 
time methods it can serve as the basic building block 
for a kinetic Monte Carlo algorithm simulating a general 
many-particle reaction-diffusion system. 

The basic idea is to generate trial moves with the ex- 
actly known single particle Green's function for a spa- 
tially constant annihilation rate k m , which is the max- 
imum of fc(r) in the current protecting domain. With 
probability k(r)/k m the particle is annihilated at the trial 
position r, otherwise a new trial move with initial posi- 
tion r is generated. The iteration proceeds until either 
the particle is annihilated or the boundary of the pro- 
tecting domain is reached. In this paper we proved rigor- 
ously the correctness of this algorithm and demonstrated 



its numerical accuracy and efficiency with an illustrative 
example. 

Important applications with a spatially varying trans- 
formation rate include continuum models for intracellular 
transport (or more generally intermittent search strate- 
gies [13|])- In intracellular transport particles (proteins, 
organelles) can switch between free diffusion and ballis- 
tic motion by molecular motor assisted movement along 
cytoskcleton filaments. The density of filaments in the 
space direction 17, pn(r, t), is generally very inhomoge- 
neous in space and sometimes even varies over time (for 
instance during cell polarization). This situation can be 
described by the Fokker-Planck equation for the proba- 
bility densities Po(r, t) and Pn(r, t) for diffusing particles 
and particles that move with a constant velocity vq in 
direction Vt, respectively 

d f 
-P M) = Z?AP M)- 7 PoM) J dn Pn (r,t) 

+i J dflPa(v,t) (45) 

J^-PfiM) = -V-(v n Pn(r,t))+7ft 1 (r,t)P (r ) t) 

-7'-Pn(r,t), (46) 

where 7 and 7' are the attachment and detach- 
ment rates (to and from filaments), respectively. The 
freely diffusing particle sees a total annihilation rate 
k(r,t) = 7/ df2 pn(r,t), with which it is transformed 
into a ballistically moving particle with a randomly 
chosen direction f2 (and velocity vq) with probability 
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pn(r,t)/ J dflpfi(r,t). The algorithm presented in this 
paper handles a Monte Carlo simulation of the diffusion 
process described by (|45|) . whereas the implementation 
of the ballistic motion (fit)]) is straightforward. 



Appendix A 

This appendix presents some analytic solutions of Eq. 
(|5|), which have mostly been taken from [l2|. Further- 
more, it derives expressions for sampling according to the 
densities p® , p® , p® in cases where this might not be ob- 
vious anymore. We list only frequently used domains in 
one, two and three dimensions. 



1. Particle on the interval [0, L] 



absorbing on both sides: 



P D (x,t\x ,t ) = j fc " 2D(t to) sm(k n x)sm(k n x ) 



with k n = r -f. 

Expressions for the probability densities p® , p® , p® 
and the corresponding distribution functions , Fj? , 
can be derived analytically. 

• reflecting on the left and absorbing on the right side: 



P D (x,t\x ,t Q ) = - kn2 ° {t ta) cos(k n x)cos(k n x ) , 

71=0 

with k n = 

Expressions for the probability densities p® , p® and 
the corresponding distribution functions F® ', F® are 
analytically derivable. As there is only x = L for the 
particle to leave the domain, it follows p®(0\t, xo, to) = 
and pf(L\t,xo,t ) = 1. 

• reflecting on both sides: 

P D (x, t\x ,t ) = 



\ 71=0 



with k n = ?f. 



* ' cos (k n x) cos (k n xo) 



Expressions for the probability density p® and the cor- 
responding distribution function F® can be derived an- 
alytically. 



2. Particle in a rectangle [0, a] x [0, b] and in a 

cuboid [0, a] x [0, b] x [0, c] 

If the boundary conditions do not vary along each side, 
Pd factorizes: 

P D = P^(x,t\x ,t ) ■ P b D (y,t\y ,t ) (-P^(z,t|z ,t )) , 



where Pg, P|> (Pf>) are given by solutions for intervals 
from the subsection above. Depending on the bound- 
ary conditions, p® is sampled by generating a random 
time for every coordinate, where there is at least one ab- 
sorbing boundary. The smallest of these times has to 
be returned as %. The particle reaches the boundary in 
the corresponding coordinate. All other quantities are 
sampled as above. 



Particle in a circle of radius R 



• absorbing boundary: 



P D (r,(p,t\r ,ipo,t ) 



ttP 2 



X cos(n(^-^ )) 



n— — oo 



j> n {a n y 



where ^ Q denotes the infinite sum over all positive 
roots a n of the Bessel function J n (a n ) = 0. 
The density of finding the particle at an arbitrary angle 
at radius r is then given by 



p r {r,t\ro,to) 



P D (r,ip,t\r ,ipo,to)rdip 



2 ag £(^ Jo(a %)Jo(a %) 

R 2 ^ J!(a ) 2 



and the corresponding distribution function is given by 
F r (r,t\r ,t ) = / dr' p r (r' ,t\r ,t ) 



2_ e _ a 2£C^oi r Ji (oof) Jo (a %) 
a Ji(a ) 2 



Hence, the distribution function belonging to p® is 
given by 



F t D (Kto) = l-2^e- a »- 



a Ji(ao) 



Analytic expressions for all quantities depending on 95 
which are needed, are straightforwardly derivable by in- 
tegrating the cos-functions. 

Having precomputcd the values of a n , random num- 
bers are sampled by inverting the occurring distribu- 
tion functions numerically. 

For a particle starting in the center of the circle the in- 
dependence becomes uniformly distributed in the inter- 
val [0, 2ir[ and F® simplifies to 

* ^ 2 O(t-tn) 

F b D (t\0,t o ) = l-2^e- a °-V L 



a Ji(a ) 



which is used to derive Eq. (fT^|) . 
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• reflecting boundary: 



P D (r,ip,t\r ,<po,t ) 



nR 2 



1+ ^2 cos (n (ip - tpo)) 

n=— oo 



Jn(^n) 2 



where ^ a denotes the infinite sum over all positive 
roots a n of J^(a n ) = 0. 

The density of finding the particle at an arbitrary angle 
at radius r is then given by 



Pr(r,t\r ,t ) 



2 

R? 



/j(r, <,j,t|ro,¥>o)rG^ 

-ag5^ni r Jo (aoj) Jo (t*o%) 
Jo(ao) 2 



and the corresponding distribution function is given by 
F r (r,t\r ,t ) = 

-al^gsl r Ji (goj) J o (<*o%) _ 
aoJo("o) 2 



r 2 2 x ^ 



A distribution function for the angle ip under the con- 
dition of being at radius r can be derived straightfor- 
wardly by integrating the cos-functions. Having pre- 
computed the values of a n , r and ip are sampled by a 
numerical inversion of the distribution functions. For 
a particle starting in the center of the circle the in- 
dependence again becomes uniformly distributed in the 
interval [0, 2tt[. 



4. Particle in a sector of angle <d with reflecting 
boundaries 



PD(r,(p,t\r 0> (po,to) 



g 4D(t-t ) 

2Q D(t - t ) 



TTip 



rr 



2D(t-t ) 



n<Po\ T frr \ 
cos n — cos n 

e / V e / * \2Dt) 



71=1 

where I u denotes the modified Bessel function of order 

tjj. 

The density for finding the particle at an arbitrary an- 
gle ip at radius r is then given by 

^2 i 2 
_ r +r O 
f g 4D(t— to) I T Tq 



Pr(r,t\r ,t ) = 



2D(t - t Q ) 



Io 



2D(t - t Q ) 



As there is no analytic expression for a distribution func- 
tion of p r available, the usage of the inversion method 
would be very slow, as the integration of r would have to 
be done numerically. Fortunately, p r does not depend on 
the sector angle 0, hence the analytically known solution 



for = 7r (half-plane) can be used to generate r for all 0. 
A distribution function for the angle ip under the condi- 
tion of being at radius r can be derived straightforwardly 
by integrating the cos-function. 



5. Particle in a sphere of radius R with absorbing 
boundary conditions 



P D (r,p,d,t\ro,po,^o,to) 
1 



2TtR 2 JT¥o~ 



£)(2n + l)P n r>(*>,tf,¥>o,0o)) 



E 



_ a 2 D(t-t a) J n+ i («n^) J n+ i 



where denotes the infinite sum over all positive zeros 
a n of the Bessel function J n+ i(a n ) = 0, P n is the n-th 
Legendre polynomial and fi is the cosine of the angle 
between r and ro- 

The density of finding the particle at arbitrary angles 
ip, at radius r is then given by 

p r (r,t\r ,t Q ) = 

dp / d$P D (r, ip, tf,t\r ,tp ,$o,to) sin(i?) r 2 
o Jo 



2rf 
R 2 ^o 



J' L (a ) 



n— 1 



2 2 D(t-t ) t /n7rr \ . /mrr\ 

^2 em I 1 Qln 1 I 



, sin 

V R J \ R J 



The corresponding distribution function can be derived 
by integrating the sin-functions: 



2 -» 2 7r 2 D(t t 0> • f n7Tr 0\ 

= — — } e sin — — 

Rr ^ V R J 

n—l 

R 2 /"nirr\ R fnirr 

— - — — sin r cos 

n 2 7r 2 V R J mi \ R 



F r (r,t\r ,t ) 



Hence, the distribution function belonging to is 
given by 



*b(t\ro,to) = 



2R^ _ n v»£fc^fll . (rmrp^ 

1 > e r 1 sm 

nr ^ V R J n 



which was used to derive Eq. (fl5|) with the help of 
l'Hospital's rule (r -> 0). 

A distribution function for /i£ [— 1; 1] under the condi- 
tion of being at radius r can be derived straightforwardly 
by integrating the Legendre polynomials P n . Using the 
sampled p,, the angels ip and $ are sampled. 
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